TO APPEAR IN ApJ 

Preprint typeset using I^T^X style cmulatcapj v. 6/22/04 



WEAK LENSING BY LARGE-SCALE STRUCTURE WITH THE FIRST RADIO SURVEY 
Tzu-Ching Chang 1,2 , Alexandre Refregier 3,4 , David J. Helfand 1 

1 Department of Astronomy and Columbia Astrophysics Laboratory, Columbia University, 550 W. 120th Street, New York, NY 10027 

2 Department of Astronomy, University of California at Berkeley, 601 Campbell Hall, Berkeley, CA 94720; tchang@astro.berkeley.edu 

3 Service d'Astrophysique, CEA/Saclay, 91191 Gif sur Yvette, France and 
4 Institute of Astronomy, Madingley Road, Cambridge CB3 OHA, UK 
To appear in ApJ 



o 
o 

(N 
bJQ' 

< 

o 

m 



> 
oo 
^t- 
in 
oo 
O 

o 
6 

c3 



ABSTRACT 

We present the first measurement of weak lensing by large-scale structure on scales of 1 — 4 degrees 
based on radio observations. We utilize the FIRST Radio Survey, a quarter-sky, 20 cm survey produced 
with the NRAO Very Large Array (VLA). The large angular scales afforded by the FIRST survey 
provide a measurement in the linear regime of the matter power spectrum, thus avoiding the necessity 
of applying uncertain non-linear corrections. Moreover, since the VLA interferometer has a well- 
known and deterministic beam, our measurement does not suffer from the irreproducible effects of 
atmospheric seeing which limit ground-based optical surveys. We use the shapelet method described 
in an earlier paper to estimate the shear from the shape of radio sources derived directly from the 
interferometric measurements in the Fourier(uu) plane. With realistic simulations we verify that 
the method yields unbiased shear estimators. We study and quantify the systematic effects which 
can produce spurious shears, analytically and with simulations, and carefully correct for them. We 
measure the shear correlation functions on angular scales of 0.5° — 40°, and compute the corresponding 
aperture mass statistics. On 1° — 4° scales, we find that the B-modes are consistent with zero, and 
detect a lensing E-mode signal significant at the 3. Oct level. After removing nearby radio sources 
with an optical counterpart, the -E-mode signal increases by 10-20%, as expected for a lensing signal 
derived from more distant sources. We use the E-mode measurement on these scales to constrain the 
mass power spectrum normalization og and the median redshift z m of the unidentified radio sources. 
We find a$(z m /2) - 6 ~ 1.0 ± 0.2 where the Ict error bars include statistical errors, cosmic variance, 
and systematics. This is consistent with earlier determinations of as from cosmic shear, the cosmic 
microwave background (CMB) and cluster abundance, and with our current knowledge of the redshift 
distribution of radio sources. Taking the prior ct§ = 0.9 ± 0.1 (68%CL) from the WMAP experiment, 
this corresponds to z m = 2.2 ±0.9 (68%CL) for radio sources without optical counterparts, consistent 
with existing models for the radio source luminosity function. Our results offer promising prospects 
for precision measurements of cosmic shear with future radio interferometers such as LOFAR and the 
SKA. 

Subject headings: cosmology: large-scale structure of universe - cosmology: dark matter - gravitational 
lensing - techniques: interferometric 



1. INTRODUCTION 

Weak gravitational lensing by large-scale structure, or 
'cosmic shear', has emerged as a powerful tool for mea- 
suring the mass distribution of the Universe (see van 
Waerbeke & Mellier 2003; Refregier 2003 for recent re- 
views). This effect is based on the distortion induced 
in the images of background galaxies by the gravita- 
tional tidal field of the intervening large-scale structure. 
Since light rays from nearby background galaxies travel 
through similar mass inhomogeneities along the way, the 
resulting image distortions are correlated. By measur- 
ing the coherent lensing distortions of background galax- 
ies, one can thus probe the mass distribution projected 
along the line of sight. The underlying physics is well- 
understood, and the observational signatures can be di- 
rectly compared to theoretical predictions. The distor- 
tion depends on the mass fluctuations, as well as on the 
mean mass density of the Universe. Since it is a projected 
effect, the amplitude is also sensitive to the geometry of 
the Universe and on the background source distribution 
in redshift space. 

The lensing distortion can be quantified by a 2 x 2 



tensor field, parameterized by the convergence and shear 
parameters, with the latter being the direct observable. 
Since galaxies have an intrinsic shape distribution, the 
shear signal must be measured statistically, assuming 
galaxies are on average randomly oriented. The lensing 
effect of interest here is not associated with any specific 
mass concentrations, but is measured by sampling ran- 
dom regions of the sky. The typical amplitude of the 
shear in this regime is only on the order of a few per- 
cent on arcminute scales, and is therefore challenging to 
measure. Large areas of the sky must be surveyed to 
reduce statistical errors, and systematic effects that spu- 
riously distort the shape of the observed sources need to 
be carefully corrected for. 

The theoretical basis for cosmic shear was pioneered 
by Gunn (1967), and was further elaborated in a mod- 
ern cosmology language by, e.g., Blandford et al. (1991), 
Miralda-Escude (1991), Kaiser (1992, 1998), Jain and 
Seljak (1997), and Schneider (1998). Because of its chal- 
lenging observational nature, cosmic shear was firmly de- 
tected for the first time only recently (Bacon et al. 2000; 
Kaiser et al. 2000; van Waerbeke, et al. 2000; Wittman 
et al. 2000). Since then, the cosmological origin and 
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implications of the shear signal have been rapidly estab- 
lished by many more groups (most recently Bacon et al 
2003; Brown et al. 2003; Hamana et al. 2003; Hoekstra 
et al. 2002; Jarvis et al. 2003; Refregier et al. 2002; van 
Waerbeke et al. 2002). Collectively, the shear two-point 
statistics have been measured from 9 ~ 1' to 100' angular 
scales, spanning the non-linear (9 < 10') and quasi-linear 
(6 > 10') regime of the matter power spectrum, and 
are consistent with predictions from the popular ACDM 
model. Constraints on the matter density Q m and the 
power spectrum normalization og from different cosmic 
shear measurements are broadly in agreement with one 
another, and are roughly consistent with measurements 
from the cluster abundance method (e.g., Pierpaoli et al. 
2001; Seljak 2002). 

Until now, all current cosmic shear measurements have 
been performed in the optical and near-infrared bands, 
with similar observational and analysis techniques. Here, 
we report a cosmic shear measurement using the FIRST 
Radio Survey (Becker et al. 1995; White et al. 1997). 
FIRST is a quarter-sky radio survey at 1.4 GHz con- 
ducted using the Very Large Array (VLA) in its B- 
configuration, and therefore provides a unique measure- 
ment of the mass power spectrum on large angular scales 
(Kamionkowski et al. 1998; Refregier et al. 1998; Chang 
& Refregier 2002). Radio surveys like FIRST offer sev- 
eral advantages for cosmic shear. Firstly, unlike optical 
galaxies, bright radio sources are typically at high red- 
shifts, thereby increasing the path-length to each through 
the Universe, and thus the strength of the lensing signal. 
Secondly, the VLA is a radio interferometer and thus 
has a well-known and deterministic beam, allowing for 
modeling of crucial systematic effects to high accuracy. 
In contrast, ground-based optical surveys are limited by 
the irreproducible effects of atmospheric seeing. Finally, 
FIRST is a sparsely-sampled but wide-angle survey, well- 
suited for measurements in the linear part of the mass 
power spectrum. This avoids the theoretical uncertain- 
ties arising from the non-linear corrections of the mass 
power spectrum required to interpret the optical cosmic 
shear surveys that are sensitive to smaller scales. 

To deal with the important issue of accurate shear 
measurements and systematics corrections, we use the 
shapelet method described in Chang & Refregier (2002, 
hereafter CR02; see also Refregier 2003b; Refregier & Ba- 
con 2003) to measure the shape and shear information 
from interferometric data directly in Fourier space. The 
method is linear and yields unbiased shear estimators. 
Systematic effects such as instrumental distortions can 
be closely modeled. 

Our paper is organized as follows. In Sj21we describe the 
theoretical background and notations for cosmic shear. 
In we summarize the main features of the FIRST 
radio survey and the properties of its radio sources. In 
HI we describe briefly our shape measurement method, 
discussing the impact of systematic effects in fJSJ The 
systematics corrections are described in Sj^J In ^3 we 
present our results for the measurement of cosmic shear, 
and discuss their cosmological implications. Our conclu- 
sions are summarized in SjS] 

2. THEORY 

We first briefly describe the theoretical basis of cosmic 
shear and the statistics we will use to present our results 



(see, e.g., Bartelmann & Schneider 2000 for a review). 

2.1. Weak Lensing 

In a weakly inhomogeneous FRW Universe, light rays 
from distant sources are deflected by the gravitational 
tidal field of intervening structures. In the weak-lensing 
limit, the resulting distortion and magnification can be 
quantified by a 2 x 2 tensor field, ^; m , which, following 
from the geodesic equation, can be related to the gravi- 
tational potential <J>: 

= ~a J d Xg(x)did m $ (1) 
where the radial window function g(x) is defined as 

9(X) = r( X ) f" d X 'n(x') r ^ { ~ ) X) . (2) 

Here, n(x') is the normalized source redshift distribution, 
r(x) is the comoving angular-diameter distance, x is the 
radial comoving coordinate, and Xh corresponds to the 
horizon. The comoving derivatives di are perpendicular 
to the line-of-sight. The symmetric distortion tensor i\)i m 
can be parameterized by the convergence, k, and the 
shear 71 and 72. These are defined as 

K = 5(^11 +V>22) (3) 

71 = ^(^11 - ^22)1 72 = ipu = "021 (4) 

Adding to Equation (0 a ^33 term which cancels out 
upon x integration, and using Poisson's equation in the 

form of V 2 $ = 3H 2^ m 5, one relates the convergence field 
k to the density fluctuation 6, weighted by the window 
function <?(%). Applying Limber's equation in Fourier- 
space, the power spectrum of convergence, P K , can then 
be expressed in terms of the 3-D mass power spectrum, 
P s (e.g., Kaiser 1992) as 

p a) - r dx \ 9{x) 1 2 c — ■ x\ 

where S is the mass-density fluctuation, a(x) is the cos- 
mological scale factor, and H and f2 m are the present- 
day Hubble constant and matter density parameter, re- 
spectively. Note that the convergence and shear fields 
are related by (e.g., Kaiser & Squires 1993) 

k = d~ 2 didjjij, (6) 

where d~ 2 is the inverse 2D Laplacian operator 

<T 2 = — / d 2 f\n\f-?\. (7) 

27T J 

In the flat sky approximation, the convergence power 
spectrum is equal to the shear power spectrum P 7 = P K . 
The shear power spectrum is simply related to other 2- 
point shear statistics which are more convenient in prac- 
tice and which we describe below. 

2.2. Shear Correlation Functions 

For each pair of galaxies, we define the tangential and 
45°-rotated shear components, 74 and j r , with respect to 
the great circle connecting the two galaxies: as 

7t=7i cos(20) + 72 sin(20) 

7 r = 7 2 cos(20) - 7isin(2(9), (8) 
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where 6 is the position angle between the x-axis and 
the great circle. In analogy to the CMB polariza- 
tion correlations, one can then construct coordinate- 
independent correlations from the rotated shear compo- 
nents (Kamionkowski et al. 1998). There are three inde- 
pendent pair- wise shear correlation functions, C\(6) = 
(7t(0o)7t(0o + 0)>, C 2 (6) = (7r(0o)7r(0o + 9)), and 
C 3 (0) = (7t(0 o )7r(0o + #))• The first two correlation 
functions are related to the shear power spectrum (e.g., 
Miralda-Escude 1991) 

c ^°) = / ^7(0Pb(W) + mw)} 

C ^0) = J l -^P 7 (l)[Jo(l9) - MW)], (9) 

while the parity invariance of weak lensing ensures that 
Cs(9) vanishes. The measurement of a non-zero Cs(0) 
is, therefore, an indication of a non-lensing contribution, 
such as that from residual systematics. 

2.3. M ap Statistics 

The aperture mass is defined as (Kaiser 1995; Schnei- 
der et al. 1998) 

M ap (9) = f £*V17(M)k(#, (10) 
J\<f>\<6 

where U{4>) is a compensated filter defined so that 

J d<j) <t>U{4>) — 0, and n{4>) is the convergence field. The 
M ap statistic thus measures the spatially-filtered pro- 
jected density field, and can be conveniently related to 
the shear as 

M a p{9) = [ d 2 4>Q{\ct>\) lt {<t>), (11) 

J\<t>\<6 

where Q{<j>) = ^ ' d$ '(j/U(<f/) ~ U ((/>), and 7t is the 
rotated tangential shear, as in Eq.ljSJ, with respect to the 
aperture center. The aperture mass variance is related 
to the shear power spectrum by (Schneider et al. 1998): 

(M 2 ap (9)) = ^ J ^P,(1)[M19)?- (12) 

The Map statistic has the advantage that the mass den- 
sity can be directly obtained from the observables - the 
shear - without the need for mass reconstruction. The 
weight function Q is relatively narrow in Fourier space, so 
that the M ap (9) measurements do not strongly correlate. 
As a result, M ap (6) and M ap (26) are almost independent 
of each other. Note that M ap {6) probes the scale of ~ 9/4 
(Schneider 1998). 

Weak gravitational lensing arises from scalar pertur- 
bations to the space-time metric, and therefore the shear 
field is expected to possess no handedness. One can de- 
compose the shear tensor field into gradient (E-mode) 
and curl (B-mode) components (Stebbins 1996). The 
gradient part contains the weak lensing signal, while the 
curl part is expected to be zero. Gravitational waves 
produce non-zero B-mode signals, but their effect is ex- 
pected to be very small. Thus, the B-mode provides a 
useful check for any non-lensing contaminations, such as 
residual systematic effects or intrinsic shape correlations 
(e.g., Heavens 2001). A rotation of 45 degrees in the 



aperture brings a curl-free field to a curl field, and also 
transforms a tangential shear to a radial shear compo- 
nent. Thus, the B-mode of the aperture mass, Mj_, is 
simply 

M±(0) = [ d 2 ^Q(H) 7 ,(0), (13) 

J\<f>\<$ 

where 7 r(</>) is the radial shear with respect to the aper- 
ture center. The aperture mass statistic therefore pro- 
vides a convenient way for E- and B-mode decomposi- 
tion, and can be computed directly from the observed 
shear. 

The aperture mass variance can also be expressed in 
terms of the shear correlation functions. Defining 

C+(0) = Ci(0) + C 2 (0); C-(0) = C 1 (9)-C 2 (9), (14) 

the aperture masses are given by 

«W> = \ [ ^1<7 + (*)T + (|) + CL(*)T_(J)] 

(Ml(9)) = \ [ ^[C + (*)T + $ - C-W-§)\ 

(15) 

where the functions T + and T_ are defined in Schneider, 
van Waerbeke & Mellier (2002). 

3. DATA 

3.1. FIRST Survey 

The FIRST radio survey (Becker et al. 1995; White 
et al. 1997) was conducted at the NRAO Very Large 
Array (VLA) at 1.4 GHz in the B configuration. The 
survey is the radio equivalent of the Sloan Digital Sky 
burvey covering ~ 10, 000 deg 2 of the Northern Galac- 
tic Cap. It consists of 3-minute snapshots covering a 
hexagonal grid using 14 3-MHz frequency channels. Its 
5cr flux density limit is 1 mJy, with a restoring beam 
FWHM of 5". 4. The survey contains ~9x 10 5 sources, 
roughly 40% of which are resolved. The basic source 
information is described in the on-line FIRST catalog 
(http://sundog.stsci.edu/). 

3.2. Redshift Distribution 

For our weak lensing measurement, the source redshift 
distribution determines the weight function along the line 
of sight, and thus conditions the strength of the lensing 
signal. The signal depends most strongly on the source 
median redshift z m and somewhat on the redshift distri- 
bution. As discussed below, the redshift distribution for 
radio sources is rather uncertain, and we have therefore 
considered z m as a free parameter throughout the paper. 

For our purposes, we consider the redshift distribu- 
tion for radio sources estimated by Dunlop & Peacock 
(1990; hereafter DP). These authors estimated the ra- 
dio luminosity functions (RLF) for the steep- and flat- 
spectrum sources at 2.7 GHz, using faint sources with 
optical counterparts in the UGC catalog and the Parkes 
selected regions database with source flux densities > 100 
mJy. They presented seven redshift models, with model 
7 considered with two cases: mean- and high-redshift 
source distributions. The former have assigned redshifts 
that are the mean values for galaxies of a given K-band 
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luminosity, and the latter have redshifts larger than the 
average, to account for a possible bias. 

The normalized differential rcdshift distribution of the 
DP models is shown in Fig. ^ The RLF's were in- 
tegrated over luminosities that produce flux densities 
greater than f mJy. The DP models were shifted to 
1.4 GHz for computation, assuming a spectral index of 
a = —0.85 for steep-spectrum sources and a = for 
flat-spectrum sources, where L a t/ a . The DP model 1 
is the fundamental model, and models 2 to 5 are varia- 
tions relative to it. For clarity, only the averaged values 
of models 2 to 5 are shown. The high redshift end was 
truncated at z = 4.5, as models 4 and 5 have unphysi- 
cal spikes at z ~ 6. Model 7 includes effects of source 
evolution in luminosity and density, and is considered 
the most probable approximation (Magliocchetti et al. 
1999). The spike near z — is due to starburst galaxies; 
however, optical follow-up observations suggest the spike 
is probably unphysical (Magliocchetti et al. 2000). The 
vertical bars show the median redshifts for each model 
(truncated at z = 4.5) which range from z m ~ 0.9 to 1.4. 
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Fig. 1. — Normalized differential redshift distribution of the 
FIRST radio sources using the models by Dunlop and Peacock 
(1990, DP). The DP model 1, the average of models 2-5, and models 
7 with the mean and high-redshift estimations are shown. The 
vertical lines indicate the median redshifts of each model. 

4. SHEAR MEASUREMENT 

4.1. Shapelet Method 

We summarize the relevant components of the shapelet 
method described in CR02. In this approach, the surface 
brightness /(x) of an object is decomposed as 

/(x) = ^/ n B n (x;/3), (16) 

n 

where the Fourier transform of the basis functions are 
given by 

5 n (x; p) = H ni (^ Xl ) H n2 (^) ^ 

[ 2 (ni+n 2 ) 7r ^2 U2 J] 2 



representing the two-dimensional orthonormal Gauss- 
Hermite basis functions of characteristic scale (3, where 
H m (£,) is the Hermite polynomial of order m, x = [x\, x%) 
and n = (ni,Ti2). In practice, the series only include a 
finite number of terms with order ri\ + n-i < iV max . 

The decomposition and its applications have also been 
studied by Refregier (2003), Refregier & Bacon (2003) 
and independently by Bernstein & Jarvis (2002). The 
bases are complete and yield fast convergence in the 
expansion if /? and x = are close to the size and 
location of the object, respectively. The basis func- 
tions can be thought of as perturbations around a two- 
dimensional Gaussian, so that the coefficients are the 
Gaussian-weighted moments of the source. The basis 
functions are also the eigenfunctions of the Quantum 
Harmonic Oscillator. Under a linear coordinate transfor- 
mation, such as a weak shear, the basis functions have 
analytic responses, making an unbiased shear estimator 
easy to construct from the coefficients. Furthermore, the 
basis functions are their own Fourier transforms up to a 
rescaling factor (see Eq. ^| below) , and are thus local- 
ized functions in both real and Fourier spaces. They are 
thus convenient for analyzing data from interferometers, 
which measure the Fourier-transform of the sky surface 
brightness. 

Making use of the basis functions, we model the source 
intensity directly in the Fourier (or uv) space, where the 
FIRST interferometric data (visibilities) are collected. 
The Fourier transform / s (k) of object intensity f s (l,m) 
of each source s is a sum of the Fourier shapelets basis 
functions 

/ s (k) =Y / fnsB n (k-k s ;fc 1 ), (18) 

11 

where _ 

B n (k;P~ 1 ) = i^ +n ^B n (k;p). (19) 

A x 2 fit is used to simultaneously model all sources in 
the field: 

X 2 = (V-M/) T CT 1 (V-M/), (20) 

where V = {V c } is the binned visibility data vector, 
M = {B c } is the theory matrix whose components are 
the basis functions of different order n evaluated at the 
respective source k s , and f = {/ ns } is the coefficient 
vector for all sources. We assume the data error matrix 
C is diagonal and constant, C = er 2 I, where a is the 
noise level and I the identity matrix; this is a reasonable 
approximation for the VLA observations since the noise 
of antennas is independent. To find the best-fit solution, 
we solve a set of simultaneous least-squares equations. 
The covariance matrix of the best-fit coefficients f is then 

W(f) = cr 2 (M T M) -1 , (21) 

which provides an estimate of the errors on the coef- 
ficients and can be used to compute the errors of de- 
rived quantities, such as source flux density and size, and 
the shear estimator. In H4.2I below, we explain how the 
shapelet centroid, scale 0, and maximum order A^ max are 
chosen in practice. 

An unbiased shear estimator for each source and for 
each shapelet order n is then given by (Refregier & Bacon 
2003) 



[n(n + 2)] 2 (f n -2,o - /n+2,o) 
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where 7 = 71 + 172 is the complex shear, and the brack- 
ets denote an average over an unlensed object ensemble. 
The coefficients f n>m are the complex polar shapelet co- 
efficients which are related to the Cartesian coefficients 
of equation i|16|) by a simple linear transformation (see 
Refregier 2003b). In this paper, we will only consider 
the n = 2 shear estimator 72, which captures most of 
the shape information of faint radio sources. The error 
on the shear estimator can be straight-forwardly propa- 
gated from the covariance matrix. 

4.2. Choice of Shapelet Parameters 

One issue which needs to be resolved is the choice of 
the shapelet parameters - (3, N max and centroid - for 
each source. The shapelets centroid was simply chosen 
to match the source centroid position in R.A. and Dec. 
listed in the FIRST catalog, converted into /, m, n coor- 
dinates. 

To generate this conversion, we note that the R.A. and 
Dec. coordinate system, also called the x, y, z coordinate 
system, is fixed to the Earth. It is defined so that the 
x, y axes define the plane of the equator and the z axis 
points along the rotation axis of the Earth. The l,m,n 
system is identical to the u,v,w system, and is defined, 
for a given pointing, so that the w axis points towards 
the phase tracking center and the u axis is in the x, y 
plane. Let h and 5 be the hour angle and declination 
of a source. Its coordinate in the x, y, z system is then 
x = cos (5 cos /i, y — — cos (5 sin ft, z = sin 8. The corre- 
sponding coordinates in the u, v, w coordinate system are 
then given by 

I \ / sinft cosft \ 

m I = I — sin So cos ho sin So sin ho cos So 
n J \ cos So cos ho — cos So sin ho sin So J 

(cos 8 cos h \ 
— cos S sin h ] , (23) 
sin S J 

where ho and So are the hour angle and declination of the 
phase tracking center for this pointing. If the source is at 
the phase tracking center (ft = ho, 8 = So), the position 
on the u, v, w system is 1 — (0,0,1), as expected since 
I, m, n are directional cosines such that I 2 + rn 2 +n 2 = 1. 
For small displacements h = ho + Ah, S = So + AS away 
from the phase tracking center, the l,m position takes 
the familiar form 1 ~ li = (— Ah cos ho, AS, 1) to first 
order in the displacement. 

The choice of f3 and N max is important to ensure that 
each source is faithfully modeled, that the shapelet series 
converges, and that our shape measurement is unbiased. 
To study the first two issues, we ran a series of simulated 
grid pointings with the observing conditions of FIRST, 
as described in CR02. The input sources are described by 
elliptical Gaussians with major axes, minor axes, and po- 
sition angles taken from the FIRST catalog. Note that 
this choice of source model is simplistic and facilitates 
the convergence of the shapelet series. This is however 
the standard model used in radio astronomy to parame- 
terize radio sources. In particular, this parameterization 
is adopted by the source fitting software used to gener- 
ate the FIRST catalog, and is therefore a natural model 
to use for our simulations. This choice may affect the 
exact values of beta and A max derived, but it will not 



bias our final shear estimator as long as these parame- 
ters are within the acceptable range (see Figure EJl. In 
CR02, we tested our shapelet reconstruction algorithm 
in detail and found that it also performs well for more 
complicated source models. 

With these simulation, we then considered a series 
of values of (3 and A max for one of the sources in the 
pointing, while keeping these parameters constant for the 
other sources (all centroids were also kept fixed). We also 
added in noise in the uv plane, and studied the behav- 
ior of sources with different SN ratios. The added noise 
level is equal to 0.15 mJy per beam in real space, which 
is the typical noise level for FIRST. The lowest SN ratio 
considered is 5, which corresponds to the FIRST detec- 
tion limit; we explored over two orders of magnitudes 
of SN ratios, covering the range appropriate for FIRST 
sources. From the reconstructed images, we find a mod- 
erately narrow range of parameters which yield sensible 
reconstructions. To quantify this range, we computed 
the x 2 difference between the input image /i n (l) and the 
fitted image /fit(l p ; (3, N max ) defined as 

x = o^W- ' (24) 

where the sum is over all A p i x pixels within a radius of 5/3 
about the centroid. For convenience, x 2 was normalized 
by the rms cr/m of the pixel distribution of f m . 

An example of the resulting dependence of x 2 on j3 and 
A max is shown in Figure |3 where the SN ratio was set to 
a high value of 500 in this case. The range of acceptable 
parameter values is reflected as the region with x 2 ~ 1 
on the figure, showing that x 2 is a good measure of the 
reconstruction. In this case, the best reconstruction is 
near ((3, A max ) = (1.5,4). In the presence of noise, the 
best reconstruction tends to favor smaller N max values 
while leaving the best (3 unchanged. 

We have repeated the above procedure for various 
sources in different simulated FIRST grids. This allowed 
us to relate the optimal /3 and A max values to the ma- 
jor and minor axes of the sources derived from the el- 
liptical Gaussian fits listed in the FIRST catalog. We 
thus derived an empirical relation between the observed 
shape parameters in the FIRST catalog (the deconvolved 
major and minor axes) and the optimal shapelets pa- 
rameters {(3 and A max ) for the source ensemble: (3 = 
^0.9 x (a/2.35) x (6/2.35), and A max = a/1.5-1, where 
a and b denote the FWHM of the major and minor axes. 
This 'recipe' was then applied to the real sources in the 
FIRST survey. 

To test the implementation of this recipe in our shear 
measurement method, we have performed further sim- 
ulations. Elliptical Gaussian sources were first gener- 
ated with source parameters following the distribution 
of the real sources in the FIRST catalog. These artifi- 
cial sources were then sheared with a constant 5% value 
for both 71 and 72. The corresponding visibilities with 
the FIRST survey settings were generated for over 1,000 
simulated pointings with ~ 30 — 40 sources per point- 
ing. We then used the above recipe to determine the 
shapelet parameters, and applied the method described 
in CR02 to calculate the shapelet coefficients for the ~ 
30,000 synthetic sources. The estimators for 71 and 72 of 
each source were then calculated according to Eq. I|22|l 
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Fig. 2. — Dependence of the shapelet reconstruction on the 
shapelet parameters /3 and N max . The color scale shows the \ 2 
difference between simulated image par ame ters and their shapelet 
reconstruction values, as defined in Eq. 1241 . Only a small range of 
/3 and 7V max values are acceptable and yield good reconstructions. 
For small j3 values, an increasing ]V max value quickly result in poor 
reconstructions. 



applied to different flux-size bins. We find that the input 
shear amplitudes are accurately recovered within the la 
statistical errors. We also performed a null simulation 
(with zero shear) and verified that the resulting shear 
measurement was consistent with zero. 

4.3. Implementation 

Given a set of FIRST flux- and phase-calibrated uv 
data, sources positioned within a radius of 64' from the 
phase center are selected from the FIRST Catalog. This 
large radius limit ensures bright sources in the primary 
beam side lobes are also included in the simultaneous fit 
and thus do not contaminate the fainter sources within 
the primary beam. Sources which are separated by less 
than 20" are merged into one source, as simulations show 
that treating close companions as one source yields the 
best shapelet fits (see also the discussion in M5.6fl . The 
shapelet parameters are determined according to the pre- 
scription described in H4.2I We then compute the source 
shapelet coefficients and their covariance matrix using a 
least-squares fit to the visibilities, where all sources in 
a given pointing are fitted simultaneously. In principle, 
some of the systematics that distort the shape of sources 
could be corrected for at this stage. However, to save 
computing time, we decided to correct for the system- 
atics afterwards, as described in the next section. The 
whole FIRST data set observed up to year 2001 were pro- 
cessed on the UK Cosmos Supercomputer. The total of 
about 50,000 individual pointings required about 9,000 
hours of Cosmos CPU time. 

As a first check, the output coefficients were exam- 
ined by computing the source fluxes and centroids de- 
rived from the shapelet coefficients (see CR02). Along 
with those fields flagged when the fit failed to converge, 



we have discarded about 3% of the processed sources, 
due to corrupted observational data, numerical conver- 
gence problems and, in a few cases, inadequate choice of 
input shapelet parameters. On average, the computed 
shapelet flux densities and centroids agree rather well 
with those in the FIRST catalog, which used an elliptical 
Gaussian fit to measure source intensities and positions. 
The shapelet flux densities and flux density errors were 
then combined to compute the signal-to-noise ratio of 
each source. To control systematics, we discarded mea- 
surements with observing parameters \HA\ > 4 hours 
and \Jl 2 + m 2 > 20', and only kept sources with an inte- 
grated flux density > lmJy and deconvolved major axis 
< 7". Excluding the measurements for point sources 
(deconvolved major axis < 2"), we thus were left with 
~ 3.6 x 10 5 usable shear estimators with associated mea- 
surement errors. 

Since the FIRST pointings partially overlap (see 
Becker ct al. 1995 for details), each source is observed 
from one to four times. We verified that the corrected 
shear estimators of a given source in different pointings 
are consistent within the errors. The shear estimators 
of a given source are then coadded using the square of 
the observed signal-to- noise ratio as a weight. The sky 
coverage in the northern cap is close to 8000 deg 2 , and 
the resulting resolved (major axis > 2" and < 7") source 
number density is about 20 sources deg -2 . 

5. SUMMARY OF SYSTEMATIC EFFECTS 

Because the weak-lensing signal is only on the order 
of 1%, systematic effects must be carefully accounted 
for, as they may otherwise introduce spurious shear 
correlations. The dominant systematics arise from the 
anisotropy of the synthesized beam, which is directly re- 
lated to the uv sampling, and thus on the geometry of 
the interferometric array. Note that, unlike the case with 
optical observations, systematic effects in our case are 
well-known and deterministic (albeit somewhat compli- 
cated) and can be calculated to high-degree of accuracy. 
In the following, we briefly describe the different sources 
of systematic effects and show that they are best viewed 
as functions of four observing parameters (see e.g., Per- 
ley et al. 1989; Taylor et al. 1999; & Thompson et al. 
1986 for a detailed description of these effects). 

5.1. The HA-DEC Effect 

Interferometric data are collected in Fourier space, or 
uv space, at a finite and discrete number of points. The 
Fourier-transform of the sky surface brightness is thus 
multiplied by the uv sampling function which determines 
the shape and size of the convolution beam (PSF) in real- 
space. The uv sampling function is deterministic and is 
a sum of delta functions centered on the positions 

u\ i / siniJ cosiJ \ f L x \ 

v ] = — I — sin 5 cos H sin S sin H cos S J I L y I , 
w J \ cos (5 cos .ff — cos 5 sin H sin S J \L Z J 

(25) 

where A is the observing wavelength, (L x , L y , L z ) are the 
VLA antenna spacings, and H and S are the source hour- 
angle and declination. Since the former two are fixed 
and known for a given observation, the effective beam 
can thus be calculated exactly from the source position 
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Fig. 3. — The artificial shear pattern induced by the HA-DEC 
effect. The length and direction of the plotted lines indicate the 
amplitude and orientation of the distortion, respectively. The line 
at the bottom-left indicates a distortion of 5%. 



Fig. 4. — The artificial shear pattern induced by the non-coplanar 
effect, for the case of (I, m) = (7' , 7') away from the phase-tracking 
center. The line at the bottom-left indicates a distortion of 45%. 



on the sky. This can be seen as the projection of the 
antenna plane onto the uv plane, which depends on the 
source position. 

With the knowledge of the uv sampling function, we 
place a point source at the phase-tracking center and 
generate simulated visibilities using the FIRST observing 
settings. The simulated visibilities are modeled using the 
Shapelet method, and the artificial shear of the simulated 
source due to systematics are calculated using the fitted 
Shapelet coefficients. We then repeat the procedure for 
various values of (H, 5). The resulting shape distortion is 
shown in Fig [3] The distortion is minimal at small hour 
angles and at declination near 34 degrees, corresponding 
to the latitude of the VLA, when sources are directly 
overhead. The induced distortion is ~ 5% at (H, 5) = 
(3.5 hr , —5°), roughly in the radial direction with respect 
to (H,S) = (0,34°). 

5.2. The Non-coplanar Effect 

The VLA is a two-dimensional array laid out across the 
surface of the Earth. Due to the curvature and rotation 
of the Earth, the visibilities are not strictly measured in a 
plane but in the three-dimensional Fourier space labeled 
as (it, v, w): 

^—27ri[ul+vm-\-w(^/l — l' 2 — m 2 — 1)] (26) 

where for simplicity we have ignored other systematic 
effects which we describe below. For short observing du- 
rations and for sources close to the phase tracking cen- 
ter, the w term is small and is often dropped. In this 
approximation the visibility function reduces to a two- 
dimensional Fourier transform 

V{u,v) = j j Vi ^™_ m2 I(l,m)e- 2 ^ l+ ™\ (27) 



which has the advantage of allowing fast computation. 
This approximation induces a distortion of the source im- 
ages that depends on the amplitude of \w(y/l — I 2 — m 2 — 
1)|, i.e., on the distance from the phase tracking center 
(I = 0, m — 0). As mentioned in the previous section, the 
amplitude of w depends on the source observing hour an- 
gle and declination. Thus, the shape distortion due to 
the two-dimensional approximation is a function of the 
source position (/, m) with respect to the phase center, as 
well as the observing hour angle and declination. In our 
Shapelet approach, we use the two-dimensional approx- 
imation in order to save computing time, and therefore 
the distortion must be corrected for. 

We quantify this distortion using simulations similar 
to that described in ^5.11 but with sources placed at 
various off-center positions in the simulated grid using 
equation l|26() . Fig. 0] shows an example of the result- 
ing distortion pattern, for which the simulated source is 
placed at (l,m)=(7', 7'). In this setting, the distortion at 
(if, 6) = (Q hr , -5°) is ~ 20%, and is roughly in the radial 
direction with respect to(H,S) = (0,34°). This effect is 
the dominant source of systematic effects (see H6.2(l . 

5.3. Bandwidth Smearing 

The bandwidth-smearing effect with interferometers is 
exactly analogous to chromatic aberration in an optical 
telescope. It is due to the fact that radio interferomet- 
ric observations are not monochromatic, but are instead 
done within a finite frequency band. As a result, the vis- 
ibility registered at a particular uv position is actually 
an average of the actual visibilities over a small inter- 
val of (it, v) positions extended in the radial direction. 
The delay tracking, which corrects for the difference in 
arrival time of light rays at two antennas prior to cor- 
relation of the signals, is only accurate at the center of 
the field-of-view and at the central observing frequency. 
Suppose that the observing central frequency is vq and 
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Fig. 5. — The artificial shear pattern induced by the bandwidth 
smearing effect. The length and direction of the plotted lines indi- 
cate the amplitude and orientation of the distortion, respectively. 



the phase tracking center is (Iq, rrio). Signals at frequency 
v arriving from sky position (l,m) will have a phase de- 
lay error of (uqI + Vom)/vQ, resulting in a phase shift of 

2tt(v — i/q)(uqI + vom)/vo- The smeared visibilities, V, 
are thus (Perley et al. 1994) 



V(u ,v ) 



dv>V[u -,v -)(-) 2 



$dv>G{v>) 

x G(v')e 2ni ^ {u ° l+V ° m \ (28) 

where G(v) is the passband function and v' = v — v a . 
For simplicity, here and in the following, we ignore 
the effect of the primary beam power pattern and the 
weighting function in the uv plane. In real space, this 
corresponds to a convolution of the true sky intensity, 
1(1, m), with a distortion function D, such as 1(1, m) — 
1(1, m) * D(l, m, v), where D is defined as 



D(l, m, v')= I / du dv e 2 * l{u ° l+v ° m) x 



dv'G(v')e 2 ^ (u ° l+v ° m) 



fdv>G(v>) 



(29) 



provided that the fractional bandwidth is sufficiently 
small. The amplitude of bandwidth-smearing distortion 
therefore depends on the source position (with respect 
to the phase tracking center), and on the size of the fre- 
quency interval. Since a change in frequency moves the 
uv points radially on the uv plane, bandwidth smearing 
distorts the observed images in the radial direction with 
respect to the field center. 

For FIRST, G(v) can be approximated by a square 
passband with a width of 3 MHz. The distortion induced 
20' away from the center of the field-of-view is about 6%. 
The resulting shear pattern is shown in Figure 

5.4. Time Averaging 



Time-averaging smearing is another systematic effect 
that results from data averaging. Correlated signals re- 
ceived by pairs of antennas are, in practice, averaged over 
small time intervals in order to reduce the size of the vis- 
ibility files. During a time-interval St, the Earth rotates 
by an angle of to e 5t, where uj e is the Earth angular ro- 
tation velocity. For observations at a declination of 90°, 
the Earth rotation corresponds to a tangential rotation in 
the uv plane. The true visibilities are therefore averaged 
over an interval of uv positions, extended tangentially, 
which results in a tangential distortion of images in real 
space. The time-averaging effect at declination 90° thus 
bears an interesting analogy to the bandwidth smearing 
effect, which distorts images in the radial direction. 

In general, however, the time averaging effect is more 
complex. Depending on the source position relative to 
the antenna (Earth) space, the loci on the uv plane over 
which visibilities are being averaged varies in length and 
shape. A convenient way to estimate the size of image 
distortion is to calculate the response of a point source. 
Since the time averaging effect preserves the integrated 
flux density, the induced tangential broadening must be 
compensated for by the reduction in the peak amplitude 
of the point source response. 

Averaging a waveform of frequency v over a time 
interval St reduces the amplitude of the response by 
smc(v5t) Ril - (m/St) 2 /6, for vdt <C 1, where v5t corre- 
sponds to the phase change. For a source at (l,m) with 
respect to the phase-tracking center, the instantaneous 
phase rate is (^l + ^ m ): an d therefore, integrating over 
a small time interval St, the reduction of the amplitude 
of the point source response is 



Rst = t 



1 - 



r 2 8t 2 



du dv 
— I + —m 
dt dt 



(30) 



where 



du uj e 

— = —(L x cos H - L y sin H) 

— ^ = —^-(L x sin 5 sin H + L y sin<5cosiJ), (31) 
dt A 

which follows directly from Equation l|25|l . and uj e = 

Combining the above two equations, we can estimate 
the reduction of response for sources at any given observ- 
ing coordinates (H, S) . The amplitude of the tangential 
width broadening due to time-averaging smearing is then 
simply Rj} . The distortion is therefore a function of the 
antenna configuration (L x and L y ), the source position in 
the sky (the hour angle and declination), and the source 
position with respect to the center of the field-of-view 
(l,m). 

For FIRST, the averaging time-interval is five seconds. 
The resulting distortion pattern for <5 = 90° in shown 
in Figure El The distortion from time averaging smear- 
ing is small in the FIRST settings; for most sources the 
distortion is less than 0.1%. 

5.5. Primary Beam Attenuation 

In practice, the true sky surface brightness 1(1, m) is 
multiplied by the antenna power pattern of the primary 
single dish A(l, m) before signals are correlated. The pri- 
mary beam power pattern is closely related to the diffrac- 
tion pattern of a circular dish, and can be approximated 
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Fig. 6. — The artificial shear pattern induced by the time- 
averaging smearing effect for the case of <5 = 90° . The distortion 
pattern is tangential in this case. The pattern varies with sky po- 
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Fig. 7. — The artificial shear pattern induced by the primary 
beam power pattern attenuation. The resulting distortion is radial 
and very small. For a big circular Gaussian source with a F WHM 
of 10", the distortion is less than 0.01%. 



by (Condon et al. 1998) 
A(r) = 



r/r 



(32) 



where J\ is the Bessel J function of order 1, r ~ 3.23 x 
Vl 2 + m 2 , and ro is the FWHM of the primary beam 
pattern. For FIRST, r = 30'.83. 

The primary beam pattern introduces a variation in 
sensitivity in the radial direction across the observing 
field-of-view which modifies the observed shape of source 
images. A general treatment of this effect is detailed in 
the Appendix. The effect depends on the primary beam 
power pattern as well as on the source sizes; for larger 
sources t he d istortion is more severe. Using Equations 
and l|A"fi|) . we calculate the effect for typical FIRST 
sources. The distortion is in the radial direction with 
respect to the field center, and the induced shear pattern 
for a circular Gaussian source with FWHM of 10" (~ 2 
times the beam size) is shown in Fig0 The distortion is 
less than 0.01% and thus the effect can be safely ignored 
in the FIRST case. 

5.6. Source Fragmentation 

A significant fraction of radio sources have a double- 
lobe structure, and are often broken into two components 
by the FIRST object finder. Because these radio lobes 
tend to be aligned, this produces strong shear correla- 
tions on small angular scales. As mentioned in ^4.31 we 
have treated sources within 20" of each other, the typi- 
cal separation of double-lobe sources, as a single source 
for the shapelet fitting. Furthermore, for sources that 
lie within 1', only one source is selected at random for 
use in the shear measurements. This leads to a loss of 
information only on scales smaller than 1', a scale much 



smaller than the scales on which we measure the weak 
lensing signal (greater than about 30'; see 33 below) 

6. CORRECTIONS FOR SYSTEMATICS 

6.1. Observations 

Without correction, the shear measurements of our 
FIRST sample are dominated by the systematic effects 
described above. These effects are complicated and are 
correlated with each other, and can not easily be decom- 
posed into convolutions and distortions, as is the case 
with optical observations. To correct them, we therefore 
adopt a statistical approach by measuring and subtract- 
ing the systematic shear as a function of the observation 
parameters (RA, DEC, 1, m) and of the source size. As 
a test, we also apply the same procedure to the simula- 
tions. 

Figs. |St and|51i show the average of the shear estima- 
tors for each of the ~ 5.3 x 10 5 measurements (with major 
axis > 2" and < 17" before coaddition) in HA-DEC and 
1-m bins for all source sizes. Since the lensing signal av- 
erages out in these cells (see M6.4I for more details), the 
resulting patterns give a measure of the systematics. 

The HA-DEC plane shows a large-scale pattern in 
15 x 6 degree bins. Most (75%) of the measurements 
concentrate in the central two HA columns, \HA\ < 1, 
which have the highest statistical accuracy. The shear 
pattern in the 1-m plane is less prominent, and shows a 
predominantly vertical pattern. The measurements are 
more evenly distributed in this plane, although the num- 
ber density decreases at larger radii. Note that all mea- 
surements with | H A | > 4hrs and Vl 2 + m 2 > 20' have 
been excluded at an earlier stage. The systematics shown 
in these two planes are coupled: the minimum shear dis- 
tortion occurs at around HA — 0, Dec ~ 34°, where 
sources are directly overhead from the VIA, and at the 
phase tracking center I = 0, m = 0. The farther away 
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Fig. 8. — (a) The shear pattern in the HA-DEC plane from the averaged shear measurements. The line at the bottom left corner indicates 
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Fig. 9. — (a) The shear pattern in the l-m plane from the averaged shear measurements. The line at the bottom left corner indicates 
an amplitude of 3% shear, (b) The shear pattern in the l-m plane from the simulated shears due to the systematic effects. The line at the 
bottom left corner indicates an amplitude of 3% shear. 



from these two central values, the greater the spurious 
shear due to systematics. We also examined the depen- 
dence of the systematics on the source size cuts, and 
found that smaller sources are subject to larger system- 
atic distortions, as expected. We do not find any sig- 
nificant dependence on observation epoch over the eight 
years of the observing campaign. 



6.2. Simulations 



As discussed in the largest systematic distortions 
for FIRST are the coupled HA-DEC and non-coplanar 
effects, which produce as much as ~ 50% artificial shear 
when a source combines low declination, with an obser- 
vation at large hour angle, and at a position far from the 
phase tracking center. Another important contribution 
is from the bandwidth-smearing effect, which produces ~ 
6% shear 20' from the phase tracking center. Although 
the systematic effects are well-understood and calculable 
to high accuracy, their impact on the shear statistics is 
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complicated, since source observations do not uniformly 
occupy the HA-DEC-l-m space and since we apply com- 
plicated cuts during the analysis. 

Ideally, we could use a control sample, such as point 
sources, to quantify the systematic effects and check our 
correction method. Since point sources carry no lens- 
ing information, any measured shear must indeed be 
due to systematic effects. However, in our shapelets ap- 
proach, the shear estimators require fitting coefficients 
up to n — 4 (Eq. 122(1 . which is a poor choice of in- 
put parameter for point sources (see, e.g., Fig. EJ. This 
not only compromises the fits for point sources, but also 
affects the reconstruction of resolved sources, since all 
sources in a given field are fitted simultaneously. Al- 
though we could imagine increasing the size of (3 for point 
sources to yield a more reasonable fit, this would be dif- 
ficult to implement given the presence of other sources 
(leaving aside the problem of quantifying the effects on 
these sources). In addition, this would have meant dou- 
bling the fitting parameters in the model and thus dou- 
bling the computing time, which was already substantial. 
More importantly, the systematic effects can in fact be 
understood and calculated very accurately, as we will 
see below. The averaged data described in the previous 
section also serves as another check for systematics. To 
quantify the systematics at the required precision, we 
therefore have carried out further simulations to study 
the impact on the shear statistics. 

For this purpose, we consider the HA-DEC, non- 
coplanar, and bandwidth-smearing effects and ignore the 
other effects, which are at least one order of magnitude 
smaller in amplitude, or have been incorporated directly 
in the analysis (such as the source fragmentation effect). 
Using the FIRST Catalog, we begin by randomizing the 
source position angles, which, along with the major and 
minor axis information, gives a measure of the intrinsic 
ellipticity distribution. For every source, we computed 



the combined systematic distortions listed above for each 
individual observation; the result is subsequently con- 
volved with the source's intrinsic shape and multiple ob- 
servations are coadded. To avoid large computing time, 
we calculated the systematics on a grid in HA-DEC-l-m 
space, and interpolate between the grid positions. Ex- 
amples are shown in Figs. |Ht> and|Hh>, where we average 
the simulated measurements in the same HA-DEC and 
1-m bins as for the observations. The resulting shear pat- 
terns and amplitudes are almost identical to those of the 
data (see Figs. |Ht and|H^ excepting for the 1-m patterns 
which show small discrepancies at the edges of the field 
where the data are of low statistical significance. We 
are therefore confident that the systematics estimation 
is, at least on average, a good representation of the true 
systematics for individual sources. 

6.3. Uncorrected Shear Correlation Functions 

We compute the rotated shear correlation functions, 
Ci, C2, and C3, as described in for both the simula- 
tion and the data. Operationally, the shear correlation 
functions are measured as 



Cx{6) = 



(33) 



where the sums are over all pairs of galaxies i and j , and 
6 = \9i — 6j\ is the separation of pairs within a bin. Sim- 
ilarly, for C3, the 7's are replaced by 7t(6 l i) and 7r(#j), 
respectively. We adopt the optimal weight function de- 
rived by Bernstein & Jarvis (2002) 

w= , 1 (34) 
N /7 2 + (l-5a J? ) 2 ' 

where o~ v is the measurement uncertainty in each shear 
component, as measured by the Shapelet method, and 
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Fig. 11. — (a) The shear pattern from the averaged shear measurements after the HA-DEC corrections. The line at the bottom left 
corner indicates an amplitude of 1% shear, (b) The shear pattern from the simulated shears after the HA-DEC corrections. The line at 
the bottom left corner indicates an amplitude of 1% shear. 



transformed to the sheared coordinates in which the 
shape is circular. We divide the data into twelve subsam- 
ples according to their sky position, and compute the cor- 
relation functions separately. The correlation functions 
are calculated in bins of Aln(0) = 0.2, and the 1-a er- 
ror bars are quantified using the field-to-field variation of 
the twelve subsamples. The error bars include noise from 
uncertainties in the shape measurements, from cosmic 
variance, and from intrinsic dispersion of the (unlenscd) 
source shapes, which by averaging over the entire sample 
is measured to be around ~ 0.4. 

The data shear correlation functions, without any cor- 
rections, are heavily dominated by the systematics, and 
we expect their correlation functions to have high ampli- 
tude and to appear similar to those of the simulation. As 
demonstrated in Figs. 110b and 110b. the general behav- 
iors of the Ci and Ci functions for the data and of the 
simulation arc in very good agreement from 10' to 40°. 
For the C3 function, the simulation captures the large- 
scale behavior well but does less well on small angular 
scales. It is nevertheless reassuring that our knowledge 
of the systematics is realistic. 

Note that the measurements of C\ and Ci on different 
angular scales in FigEJ appear correlated both for the 
data and for the simulations. This reflects the fact that 
the systematics, which dominate C\ and C2 before the 
corrections, are strongly correlated within a patch of the 
sky used to measure the correlation functions. On the 
other hand, C3 (for 9 < 1000') is less correlated, and 
therefore its errors are dominated by statistical uncer- 
tainties. 

6.4. Corrected Shear Correlation Functions 

First, let us re-examine the averaged shear pattern 
on the HA-DEC plane (Fig. Efe,). It is evident that 
this large-scale pattern is contributed almost entirely 
by systematics. Using the averaged data, we fit a two- 



dimensional 2nd-order polynomial to this plane, and sub- 
sequently subtract the fitted systematics from individual 
sources. The sources are grouped according to their sizes 
and observed 1-m position in the field for the modeling 
and the corrections. We have tried both polynomial fits 
and spline fits to the systematics and find that the result- 
ing residuals are comparable. The residual shear pattern 
is shown in Fig. 1111 and appears small in amplitude and 
random in orientation, as desired. We discard measure- 
ments with averaged residual shear values greater than 
1.5% in all following analyses. To test for our control 
of the systematics, we perform the same procedure with 
the simulated shears; the residual shear pattern is ran- 
dom and very small in amplitude. 

We then study the shear correlations function after the 
HA-DEC corrections for both the data and the simula- 
tion. As shown in Fig. 112b . the amplitudes of the data 
C\ and C2 correlations are an order-of-magnitude smaller 
than those of the uncorrected data (Fig. HOb). and the 
high amplitude tail of C3 on large angular scales seen 
before the correction has disappeared. The shear corre- 
lation functions of the simulation are also well-behaved 
on large angular scales, but have much smaller ampli- 
tudes between 10 and 100' (see Fig. EJd). 

The correction for the l-m dependence is more com- 
plicated. After source coaddition, which depends on the 
declination of the phase center and the number of times 
a source was observed, for most sources the l-m shear 
pattern roughly cancels out, while for others it is en- 
hanced or distorted. This is potentially a contamination 
to the lensing signal that is difficult to isolate and cor- 
rect for. We nevertheless perform a correction on the 
l-m plane similar to that for the HA-DEC plane. First, 
we examine the shear pattern on the l-m plane using ob- 
served and simulated shears that have been corrected for 
the HA-DEC pattern, fit a low-order polynomial to the 
plane, and then subtract the fitted pattern from individ- 
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Fig. 12. — (a) The shear correlation function of the data after HA-DEC corrections. The lines over-plotted indicate predictions from the 
fim = 0.3 ACDM model, with erg = 0.9, and source redshifts lying on the plane of z m = 1, 1.5, 2,2.5, respectively, (b) The shear correlation 
functions of the simulation after HA-DEC corrections. The corrections follow the exact same procedure as the corrections to the data. The 
lines show the ACDM predictions as in (a). 



ual measurements. The sources are grouped according 
to their size for modeling and corrections, and measure- 
ments with averaged residual shear values greater than 
1.5% are discarded. After corrections and coaddition, 
we again compute the correlation functions for both the 
data and the simulation. Interestingly, the correlation 
functions are almost identical to those that have been 
corrected for HA-DEC patterns only; i.e., the l-m pat- 
tern correction does not alter the correlation functions of 
the data or of the simulation. We also checked our cor- 
rections by considering subsamples with different source 
sizes and did not find significant differences in the result- 
ing correlation functions. 

We also ran a further series of correlation function 
tests. We consider using only sources within 15' of 
the phase-tracking center where the l-m distortions are 
small, excluding sources at low declination (< 10°) where 
the combined HA-DEC and I — m distortions are se- 
vere, relaxing or strengthening our imposed constraints 
on excluding source close companions (discarding sources 
within 5" to 60" of one another) , and changing the range 
in major axes of the selected sources. Remarkably, the 
shear correlation functions of the data and of the simu- 
lation remain very similar to those of Figs. 112b and 1 12b. 
respectively, in every case. 

We therefore conclude that the HA-DEC corrected 
shear correlation functions are fairly robust, and that 
the systematics dependence on the four parameters, HA, 
DEC, I and m are well-reproduced by the simulations. 
We thus take the simulated shear correlation functions 
after HA-DEC corrections, which are non-zero on most 
scales, as our best estimation for the residual systemat- 
ics in the HA-DEC and l-m parameter space. To correct 
for these, we subtract the residual simulated C\ and Ci 
correlation functions from the corresponding data corre- 
lation function. The errors in each correlation function 



are then added in quadrature. 

7. RESULTS 
7.1. M ap Statistics 

Using the shear correlation functions C\{9) and 6*2(0) 
measured above, we compute the M ap 2-point statistics 
using Eq. lfT31) . The variances (M| p ) and (Mf) contain 
the E-mode and B-mode, respectively, and are plotted in 
Fig. El as a function of aperture radius 9. While the E- 
mode displays a significant signal, the B-mode amplitude 
is consistent with on all scales. The scale dependence of 
the E-modes for 9 > 200' is consistent with that expected 
for the ACDM model (see curves) , but deviates from this 
model on smaller scales. The M ap statistics with an aper- 
ture radius of 9 is mostly sensitive to angular scales of 
about 0/4. The depleted E-modes may therefore reflect 
the presence of systematics on scales smaller than 50', 
which roughly corresponds to the scale of the individual 
fields (40' in diameter). This could be due to an over- 
or under-correction of one or several of the systematic 
effects discussed in section H6.4I 

In addition to the checks described in M6.4I we divided 
the data into high- and low-latitude halves, and east and 
west hemispheres, and computed again the M ap statistics 
in each case; the results are consistent with that of Fig. 
1131 the lensing signals are not confined to any particu- 
lar part of the sky. Similarly, we separated the sources 
into two groups according to the observed flux density; 
again, the M ap statistics do not differ significantly. As a 
comparison, the M ap statistics computed from the sim- 
ulations are shown in Fig. El The E- and B-modes of 
the simulation are consistent with zero on most scales, as 
expected, with the exception of slightly negative E-mode 
values on scales smaller than ~ 200' (50' in real space). 
This suggests a possible contamination from systematics 
on small angular scales, as discussed above. 
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Fig. 13. — Map statistics as a function of aperture radius 9. The 
top panel shows the E-mode variance (M^ p ), while the bottom 
panel shows the B-mode variance (M^). The curves indicate pre- 
dictions from the Q m = 0.3 ACDM model, with ct 8 = 0.9, T = 0.21, 
and several source median redshifts z m . 



In order to test for the presence of a lensing signal on 
scales larger than 50', we removed FIRST sources for 
which an optical counterpart could be identified in the 
APM catalog (McMahon & Irwin 1992), which amounts 
to 10% of the FIRST sample. Excluding quasars, which 
are both rare and mostly point sources (all of which are 
excluded from our source sample), the redshifts of the 
APM sources peak at around z ~ 0.15 and are for the 
most part bounded by z < 0.3 (McMahon et al. 2002). 
By excluding them, we increase the median redshift of 
the source sample and thus expect the lensing signal to 
increase. Fig. 1151 shows the M ap statistics for this new 
sample. The E-mode signal has a significance of 3.6a, as 
measured at the 6 ~ 450' scale. Compared to Fig. ^3 the 
E-mode signal on scales 300' < 9 < 700' is indeed larger 
by 10-20%. This confirms the presence of the lensing 
signal on these scales. 

As an interesting exercise, we calculate the median red- 
shift derived from the DP redshift models (see M3.2|) by 
excluding the z < 0.3 region; the various median red- 
shifts shown in Fig. increase by about 10-15%. Since 
the M ap lensing signal increases roughly as z^ to first or- 
der, the changes in the measured lensing signal wrought 
by the exclusion of the low-redshift sources is consistent 
with that expected from the consequent change in the 
estimated z m from the models. 

7.2. Cosmological Implications 

Using the M ap statistics from the sample without opti- 
cal counterparts, we fit cosmological models to our data 
by computing the x 2 values: 

X 2 = (d-m) T W- 1 (d-m), (35) 

where d is the M ap data vector, m is the ACDM model 
vector and W the covariance matrix. For d we use the 
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Fig. 14. — Same as the previous figure, but the M ap statistics 
are computed using the simulations, which serve as a null test. 



6x10-6 
4xl0- 6 
2x10-6 







-2x10- 



= 2.5 

= 2 

= 1.5 

= 1 



1 00 



1000 



8 [arcmin] 



6xl0- 6 
4x10-6 
2x10-6 







-2x10- 



8 [arcmin] 



Fig. 15. — Same as Fig. 1131 but this time for only those sources 
lacking optical counterparts in the APM catalog. 



M ap values with 9 > 200', thus avoiding sub-pointing 
scales which have unreliable systematics correction (see 
discussion above). The covariance matrix is computed 
from the field-to-field variations as described in t|6.3l and 
is given by 



Z n (d in - (dj)){d jn - (dj)) 
N(N-l) 



(36) 



where the sum is over different subsamples n, and N = 12 
is the total number of subsamples. The d^s are the M ap 
E-mode values at different scales, and (di) is their av- 
erage over all 12 fields. Since the B-modes quantify the 
contamination from systematics and are consistent with 
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zero at the selected angular scales, we add the B-mode 
covariance to the E-mode covariance matrix. As a result, 
the covariance matrix accounts for all the sources of er- 
rors, namely, intrinsic source shapes, shape measurement 
errors, cosmic variance, and systematics. 

The data points from angular scales 200' < 9 < 1000' 
are partially degenerate and contain only 3-4 indepen- 
dent measurements. We therefore use a singular value 
decomposition to calculate W _1 , and discard the sin- 
gular values of W which are negligible compared to the 
rest. Since the M ap data vector was initially computed 
at small angular intervals, we only keep four data points 
which contain the highest signal-to-noise ratios and are 
independent. 

The redshift distribution of our sample is rather un- 
certain, as indicated by the various models for the ra- 
dio source redshift-distribution discussed in 43. 2\ We 
therefore chose to vary two parameters in our fit: as, 
the mass power spectrum normalization in spheres of 
8/i _1 Mpc, and z m , the median source redshift. We as- 
sumed a ACDM model with fixed values of tt m = 0.3 
and r = 0.21. Note that, although we are probing scales 
greater than 8/i _1 Mpc, the parameter as is convenient 
for comparing our results with those from other groups 
and methods. 

The resulting \ 2 contour plot is shown in Fig. 1161 The 
solid contours indicate the 68.3%, 95.4% confidence levels 
from FIRST, excluding the predominantly low-redshift 
objects with APM optical counterparts. For comparison, 
the dashed contours show the 68.3% CL constraint from 
the FIRST sources including those with the APM coun- 
terparts which, as expected, are consistent with lower z m 
values. As a check, we used the linear power spectrum 
and the fitting formula from Smith et al. (2003) and 
Peacock & Dodds (1996), and found the contours do not 
vary appreciably. This is expected since we are probing 
the linear part of the mass power spectrum. To a good 
approximation, the contours (excluding APM counter- 
parts) correspond to 

a 8 (^)°' 6 -0.95 ± 0.22, (37) 

where the 68%CL error includes statistical errors, cosmic 
variance, and systematic effects. 

Recent cosmic shear measurements in the optical band 
yield values of as between 0.7 and 1.0 for O m ~ 0.3 
and T ~ 0.21 (Bacon et al. 2003; Brown et al. 2003; 
Hamana et al. 2003; Hoekstra et al. 2002; Jarvis et al. 
2002; Massey et al. 2003; Refregier, Rhodes, & Groth 
2002; Rhodes, Refregier & Groth 2004; van Waerbeke et 
al. 2002). The averaged constraint from several of these 
surveys is as — 0.83 ± 0.04 (as compiled by Refregier 
2003) for the same values of f2 m and T. The constraint, 
a s = 0.9 ± 0.1 (68%CL), from the WMAP CMB experi- 
ment (Spergel et al. 2003) is also shown in Fig. El For 
source redshifts in the range 1.4 < z m < 3.4, our results 
are consistent (within la) with both of these different 
measurements. Reversing the argument and taking the 
WMAP determination of as as a prior, we find that the 
median redshift of the radio sources in our sample (with- 
out APM counterpart) is z m = 2.2 ± 0.9 at 68%CL. This 
redshift range is also consistent with the models for the 
radio source redshift distribution described in ^13.21 

8. CONCLUSIONS 
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Fig. 16. — Constraints on the spectrum normalization, as, 
and on the source median redshift, z m , from our FIRST cosmic 
shear measurement. The solid lines indicate the 68% and 95% 
CL from FIRST, excluding the predominantly low-redshift objects 
with APM optical counterparts. The dashed lines are the 68% 
CL from FIRST including sources with the APM counterparts. A 
ACDM model with Q m = 0.3 and T = 0.21 was assumed. Also 
shown for comparison is the constraint from the WMAP CMB 
experiment <r 8 = 0.9 ± 0.1 (68%CL; Spergel et al. 2003). 



We have presented our cosmic shear measurement us- 
ing the FIRST Radio Survey over 8,000 square degrees 
of sky. We apply the shear measurement method de- 
scribed in Chang & Refregier (2002) to measure shear 
directly in Fourier space, where interferometric data are 
collected. In this shapelets approach, we carefully deter- 
mined the input parameters for the source shape decom- 
position, and verified that the shear estimators are un- 
biased and robust using realistic simulations. The dom- 
inant systematic effects associated with interferometric 
observations were studied in detail. The analytical and 
simulated expectations of the systematics compared well 
with the data when considered functions of four obser- 
vational parameters. We corrected the systematics both 
in this parameter space and in the correlation function 
domain, and verified the accuracy of the corrections with 
further simulations. 

Using the corrected shear correlation functions, we 
computed the aperture mass statistic, which decomposes 
the signal into E- and B-modes, containing the lensing 
signal and the non-lensing contributions, respectively. 
We find that the B-modes are consistent with on all 
scales considered, while the E-modes displays a signifi- 
cant lensing signal. The E-mode signal has a significance 
of 3.0a, as measured at the 9 ~ 450' scale. The signal 
increases by 10-20% when sources with optical counter- 
parts (and, therefore, predominantly at low redshift) are 
excluded. This confirms the presence of a lensing signal 
which is expected to increase as the mean source redshift 
increases. 

Using the E-mode M ap statistics on scales 200' < 9 < 
1000' (corresponding to physical scales from about 1° 
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to 4°) , we constrained jointly the power spectrum nor- 
malization as and the median redshift z m of our radio 
source sample. We find cr 8 (z m /2) - 6 ~ 1.0 ± 0.2 where 
the ler error bars include statistical errors, cosmic vari- 
ance and systematics. A ACDM model with Cl m = 0.3 
and T = 0.21 was assumed. This result is consistent 
with earlier determinations of as from cosmic shear, and 
with the WMAP CMB experiment and cluster abun- 
dance tests. Taking the prior on er 8 from the WMAP 
experiment, this corresponds to z m — 2.2 ±0.9 (68%CL) 
for radio sources without optical counterparts, consis- 
tent with existing models for the radio source luminosity 
function. 

Our shear measurement is complementary to the ear- 
lier cosmic shear results which were all performed in the 
optical or near-IR band. It is the first measurement us- 
ing radio sources observed with interferometers and is 
therefore subject to very different systematic effects. We 
used a new shear measurement method which allows di- 
rect computation of shear estimators in Fourier space. 
Moreover, our measurement probes large angular scales, 
which are in the linear regime of structure growth, elim- 
inating the need for uncertain non-linear corrections in 
the matter power spectrum. Our measurement thus of- 
fers promising prospects for measuring cosmic shear with 



future radio interferometers such as LOFAR and SKA, 
which offer exciting opportunities for precision measure- 
ments of cosmological parameters (Schneider 1999). 
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APPENDIX 

EFFECT OF SENSITIVITY VARIATIONS ON SOURCE SHAPES 

In general, the sensitivity of the image is not exactly constant across the image. Here, we compute the effect of the 
sensitivity variations on the shapes of sources. We then apply our results to the case of a radial primary beam pattern 
affecting a circular Gaussian source. 

Let z(x) be the intrinsic intensity of a source as a function of position x on the image. The observed intensity is 
i'(x) = i(x)s(x), where s(x) is the sensitivity function. In the application below, s(x) is simply the primary beam 
function. The true centroid x° of the source is defined by 



c° = J d 2 x Xii(x) J J d 2 x i(x) . 



(Al) 



We will assume that the source angular size is small compared to the scale on which s(x) varies. By expanding s(x) 
in a Taylor series about x°, it is easy to show that the observed centroid position x 0/ is, to lowest order, 

x 0/ ~x° + JijSj, (A2) 

where the summation convention was used and Jij is the normalized quadrupole moment of the source defined as 
Jij = J d 2 x (xi — x®)(xj — x { j)i(x)/ J d 2 x i(x). We have also defined the normalized derivative tensors of the sensitivity 
function as 



Si 



ds(x°) 

dxi 



A(x°) , 



Sij 



g> 2 s(x°) 

dxidx-j 



A(x°) 



We can also compute the observed normalized quadrupole moments of the source 



J'ij EE J d 2 X ( Xi - xf)( X j ~ X°'K(X) I J d 2 X i'(x) 



about the observed centroid position x . To lowest order, we find 



Jij — Jij + SkJijk ~ -^SklJklJij + ^SklJijkl ~ SkSlJjkJjl, 



(A3) 



(A4) 



(A5) 



where Jij,Jijk Jijkl are the true normalized second, third and fourth order moments of the source about x . The 
terms in Ski arise from the second order derivative of the sensitivity function which distorts the image. The terms in 
Sk arise from the shift from the true centroid. 
We will now focus on the case of a radial primary beam pattern for which the sensitivity s(x) = S(x) is only a 

function of the radius x from the center of the image. In this case, Si — ^-&», and Sij — ^-XiXj + -§^{Sij — XiXj), 
where Xi = Xi/x is the unit radial vector, and s' and s" are derivatives of s with respect to x. For simplicity we 

consider a Gaussian source of rms radius r, such that i(x) oc e" 2 ^ 1 . In this case, the true moments are simply given 
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by Jv 



J, 



tjk 



0, and Jim = J 2 



3r 4 , J u 



expressions to equation I|A5J| . we can find the observed ellipticity e£ 
<- 2 ' = \{J' lx + J'22) of the source. We find 



r 4 and Jm 2 = J12 



0. By applying these specific 



{J'n — J 22 , 2Ji 2 }/(Jii + J22) an d square radius 



r 

2s 



and 



r 

2s 



s 
x 



s H 



>\2 



( S ') 2 



(A6) 



(A7) 



where e[ = {x\ — x\, 2x1X2}/ {x\ + x^) is the unit radial ellipticity vector. Thus, the observed ellipticity pattern will 



be radial (tangential) if the quantity 
the specific case of the VLA primary beam. 



(s'f 



is positive (negative). In section M5.5I we apply these results to 
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